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Abstract 

We developed a fast numerical method for complex symmetric shifted 
linear systems, which is motivated by the quantum-mechanical (electronic- 
structure) theory in nanoscale materials. The method is named shifted 
Conjugate Orthogonal Conjugate Gradient (shifted COCG) method. The 
formulation is given and several numerical aspects are discussed. 

1 Introduction 

The quantum-mechanical feature of electrons plays a crucial role in nanoscale 
materials and its mathematical foundation is reduced to linear-algebraic prob- 
lems with given large matrices H , called Hamiltonian. The physical properties 
of electrons can be described by the Green's function G that is defined as in- 
versed matrix G{z) := {zl — H)^^ with a complex variable z whose real part 
corresponds to energy. Since the standard matrix-inversion procedure requires 
an impractical computational cost in case of large matrices or in nanoscale ma- 
terials, there is a strong need for the fast solution of the Green's function, see 
the references in 

Here we introduce a new method for calculating the Green's function [3]. 
Let H be an N-hy-N real symmetric Hamiltonian matrix, then any element of 
the Green's function can be written as 

G,,(z) = ef(z/-i/)-ie„ (1) 

where denotes the zth unit vector, and the complex energy z = a + iS E C . 
Note that the i,j entry of can be obtained by two steps: first, compute 
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{zI—H)x = Cj, and then compute ejx. Since an integral of Gy (z) with respect 
to (J is needed to obtain any physical quantity (HI, the numerical integration 
requires {(o-fc + i5)I — H^x^''^ = for k = 1,...,A/. Hence, the problem 
reduces to solving the following shifted linear systems with complex symmetric 
matrices: 

{A + akT)x^^'^ = b for k = l,...,M. (2) 

The paper is organized as follows: in the next section, we describe the algorithm 
and the property of COCG for solving complex symmetric linear systems. In §3, 
to solve efficiently, a numerical method named shifted COCG is proposed 
and seed switching technique is introduced. In §4, we report some numerical 
experiments. Finally, we make some concluding remarks in §5. 



2 The COCG method 

Matrix A is called complex symmetric if A is not Hermitian but symmetric 
A = 7^ A^ . To solve the linear systems, the COCG method 4» has been 
proposed and is known as one of the most successful Krylov subspace methods. 

Algorithm 1: COCG 

xq is an initial guess, 

ra = b- Axq, = 0, = 0, 

for n = 0, 1, . . . until ||r„|| < ei||6|| do: 

Pn = rn + f3n-lPn^l, 
T 

PlAp^ ' 

T 

r, _ '"'n+l'^'n+l 

end 

Observing Algorithm 1, we see that the nth residual can be written as r„(:= 
b — Axn) = Rn{A)ro, where Rq{X) = 1, i?i(A) = (1 — aoX)Ro{X), and 

RnW = f 1 + - a„_iA)i?„_i(A) - ^^^^Q;„_ii?„_2(A). (3) 

It is known that if breakdown does not occur, then the nth residual satisfies 

r„ ± Kn{A, ro), (4) 
which leads to conjugate orthogonality 1. rj for i ^ j. 
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3 A shifted COCG method 

In this section, we develop the COCG method for solving complex symmetric 
shifted linear systems. First, let us describe the the following theorem: 

Theorem 1 (Prommer [2, Theorem 1]) Let Wi C W2 ^ ■ ■ ■ Wk be a 

sequence of nested subspaces of such that Wn has dimension n and ]¥„ H 
(X„+i(A, b))^ = {0}, 71 = 1, . . . , fc. Let r„ := i?„(A)6, := R^{A + al)b be 
residual vectors satisfying 

r„,< _L Wn, n = l,...,k. (5) 

Then r„ and rj^ are collinear. 

Corollary 1 Let r„ and be the residual vectors of COCG started with Xq = 
a^Q = 0. Then, Vn and are collinear, i.e., there exists ttJ^ G C such that 

Proof. Since it follows ft'om 10} that the COCG residuals satisfy lO with Wn = 
Kn{A, b) = Kn{A + aL, b), this result follows from Theorem 1. 

Next, we give the formulas for computing rn+i by using the information of 
Tn+i- It follows from the polynomial JSJ that we have 

r„+i = IH Q!„ - a„A r„ a„r„_i, (6) 

<+i = + + (7) 
Substituting the relation r„ ~ '^n'^n ^^to the previous recurrence 10, we have 
[^ + ^^n-(^niA + aL) -—r„--—— r„_i. (8) 
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To obtain the computational formula for rj^^^^, the three parameters a^, Pn-i^ 
and 'Kn+i ^^"S essentially required. Hence, we give the formulas for the three 
parameters. First, comparing the coefficients of Ar„ in ^ and ©, we find 

<=(</<+i)«». (9) 

Second, comparing the coefficients of r„_i leads to ^-2— ia„ — Sub- 



stituting the result of 10 into the previous equation, we have 

(<-l/<)Vl. (10) 

Finally, comparing the coefficients of r„, we find (1 + J^^^ — a^cr) = 
1 + ^""^ a„. Substituting ||3 and (|10|) into the previous equation, we obtain 

K+l = f 1 + ^^^Cln + an<j)TTn - ^^^^a„7r,^_i = i?„+i(-(T). (11) 

V an-1 ' a„_i 
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The above formulation is based on the assumption that the seed and shifted 
systems are Ax = b and {A + al)x — b. Similarly, it can be readily general- 
ized to solve m shifted linear systems {A + (Jil)x^^^ — b using the seed system 
{A + asl)x = b. The resulting algorithm is given below. 



Algorithm 2: Shifted COCG 



■-0 



0, TT- 

for n = 0, 1 



until 



1 - a-i 
"r 



< ei do: 



Pn = +/3„-iP„_i, 



" pliA + aJ)p„' 

{Begin shifted system} 
for s) = 1, . . . , TO do 
if Wr^n^W > e2\\b\\ then 



'n+1 



is,^) 



end if 
end 

{End shifted system} 

rn+i = f„ - a„(A + (^sI)Pn, 
T 



end 



Shifted COCG with seed switching technique 

We can see from Algorithm 2 that if iTri"'*^] = \Rn\(Js — > 1, then ||ri*^|| < 
||r„||. Hence, if we could find a seed system such that \Rn\<Js — ^01 ^ 1j then 
all shifted systems could be solved. However, it is extremely hard to find such 
system in a priori except some special cases discussed in [2]. In this section, we 
will avoid such problem by using the following strategy: 

I. Choose a seed system, and then start Algorithm 2; 

H. If the seed system was solved at nth iteration, then find the new one; 

HI. Start Algorithm 2 from [n + l)th iteration using the new seed system. 

In II, as one of criteria for choosing the new seed system s, we adopt s = 
argmaxig7{||rl*'' II}, where / denotes an index set of unsolved systems. In HI, 
we need two steps to switch the old seed system to the new one. First, compute 

for obtaining r^^\i and (3n''pn\ Since it follows from rl^]_i + Pn'^pii^ that we 
have p^+i, we can start COCG solving the system {A + asl)x^^^ = b from 
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1st seed: k=301 
2nd seed: k=494 
3rd seed: k=712 ■ 



400 600 
Linear systems (k) 





1 St seed; k=501 
2nd seed: k=712 









400 600 800 
Linear systems (k) 



Figure 1: The true relative residual 2-norm histories after each seed switching 
finished. The first seed is k = 30f on the left and k — 50f on the right. 



(n + f )th iteration step. Second, to solve remaining systems by using the new 
seed s, it requires generating all\_i, pll^ from the new seed. We see that they 
can be readily generated by the following polynomial: 

■^n+i = Rn+i{(^~s - (^i) for aU i e /. 
To obtain the above polynomial, we need to compute 

for I = 0, . . . , n, J = 0, . . . , n — 1. Hence, the switching strategy requires only 
scalar operations, and moreover we can see that if breakdown does not occur, 
iterating the process from (II) to (III) enables us to keep solving the systems 
without losing the dimension of the Krylov subspace that has been generated 
until the last switching. 



4 Numerical examples 

In this section, we report the results of numerical experiments. The problem 
originally comes from 3 and is written as follows: 

{aki - H)x^''^ = ei, k = l,...,m, 

where, ak = 0.4 + (fc - 1 + i)/1000, H e ^^2048x2048 jg symmetric matrix, 
ei = (1, 0, . . . , 0)^, and m — 1001. Since {aki ^ H) is complex symmetric, 
the iterative solvers we used are COCG and shifted COCG. We can also apply 
shifted Bi-CGSTAB(^)^2_ and GMRES |lj to the above problem since they can 
be used for general non-Hermitian shifted linear systems. However, they do 
not exploit the property of complex symmetric matrix. This leads to high 
computational costs per iteration step. 

All experiments were performed on an ALPHA work station with a 750MHz 
processor using double precision arithmetic. Code were written in Fortran 77 
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Table 1: Numerical results of COCG and shifted COCG. 



Initial seed 


Switching 


Total MVs 


k = 301 


2 


330 


k = 501 


1 


328 


COCG 




124396 



and compiled with the optimization option -OA. The switching criterion is 
s — argmaxig/{||rl'^ II}. The stopping criteria are ei,e2 < 10~^^. We report 
two examples for k = 301 and k = 501 as a initial seed system. True residual 
2-norm histories are given in Fig. 1. 

In Fig. 1 on the left, 518 systems remained unsolved when the first seed 
system k — 301 converged. Then, the next seed k ~ 494 solved 412 more 
systems. Finally, the third seed k = 712 solved all of the remaining systems. 
In Fig. 1 on the right, we chose k — 501 as an initial seed system. This choice 
led to 147 unsolved systems. Finally, the next seed k — 712 solved all of the 
remaining systems. 

Numerical results are shown in Table 1. Total MVs denotes the total number 
of matrix-vector multiplications. We can see from Table 1 that shifted COCG 
required only about 0.27% of Total MVs of COCG. 

5 Concluding remarks 

Since the original problem |^ is a fundamental quantum-mechanical equation, 
the present method is applicable, in principle, to various nanoscale materials, 
such as silicon, carbon, metals, polymers and so on, if the Hamiltonian for 
electrons is given as an explicit matrix H. 

The present paper gives an interdisciplinary research between mathematics 
and physics, which shows that the computational science can give an impor- 
tant contribution to nanoscience through the development of general numerical 
algorithms, when a fundamental equation is formulated in physics. 
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